//
// ********************************************************************
// * License and Disclaimer                                           *
// *                                                                  *
// * The  Geant4 software  is  copyright of the Copyright Holders  of *
// * the Geant4 Collaboration.  It is provided  under  the terms  and *
// * conditions of the Geant4 Software License,  included in the file *
// * LICENSE and available at  http://cern.ch/geant4/license .  These *
// * include a list of copyright holders.                             *
// *                                                                  *
// * Neither the authors of this software system, nor their employing *
// * institutes,nor the agencies providing financial support for this *
// * work  make  any representation or  warranty, express or implied, *
// * regarding  this  software system or assume any liability for its *
// * use.  Please see the license in the file  LICENSE  and URL above *
// * for the full disclaimer and the limitation of liability.         *
// *                                                                  *
// * This  code  implementation is the result of  the  scientific and *
// * technical work of the GEANT4 collaboration.                      *
// * By using,  copying,  modifying or  distributing the software (or *
// * any work based  on the software)  you  agree  to acknowledge its *
// * use  in  resulting  scientific  publications,  and indicate your *
// * acceptance of all terms of the Geant4 Software license.          *
// ********************************************************************
//
//
// --------------------------------------------------------------
//   GEANT 4 - Underground Dark Matter Detector Advanced Example
//
//      For information related to this code contact: Alex Howard
//      e-mail: alexander.howard@cern.ch
// --------------------------------------------------------------
// Comments
//
//                  Underground Advanced
//               by A. Howard and H. Araujo 
//                    (27th November 2001)
//
//
// PrimaryGeneratorAction program
// --------------------------------------------------------------

#include "DMXPrimaryGeneratorAction.hh"
#include "G4GeneralParticleSource.hh"
#include "DMXParticleSourceMessenger.hh"
#include "G4GeneralParticleSource.hh"
#include "G4ParticleGun.hh"
#include "G4Event.hh"
#include "G4Run.hh"
#include "G4UnitsTable.hh"
#include "G4ParticleTable.hh"
#include "Randomize.hh"
#include "G4RandomDirection.hh"
#include "G4SPSAngDistribution.hh"

#include "globals.hh"

const G4double PI = 3.14159265358979;

DMXPrimaryGeneratorAction::DMXPrimaryGeneratorAction(G4int RunID) 
  : G4VUserPrimaryGeneratorAction(),particleGun(0)
{
  runid = RunID;
  G4int Nparticle = 1;  
  particleGun = new G4ParticleGun(Nparticle);
  //  particleGun = new G4GeneralParticleSource();
}
DMXPrimaryGeneratorAction::DMXPrimaryGeneratorAction() 
  : G4VUserPrimaryGeneratorAction(),particleGun(0)
{
  runid = 0 ;
  G4int Nparticle = 1;  
  particleGun = new G4ParticleGun(Nparticle);
  //  particleGun = new G4GeneralParticleSource();
}


DMXPrimaryGeneratorAction::~DMXPrimaryGeneratorAction() {
  delete particleGun;
  delete particle;
  delete particleTable;
}


void DMXPrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent) {

  G4int eventid = anEvent->GetEventID();
  G4int iz   = runid / 305;

  G4int nn;
  G4int NOP = 1;
  //  G4int NOP = 36;
  //  G4int NOP = 1000;

  for(int i=0;i<NOP;i++){
  //--- Ion
  //SetParticleIon( eventid );
  //--- Alpha
  SetParticleAlpha( eventid );
  //--- Photon
  //SetParticlePhoton(runid,i);
  //--- Neutron
  //SetParticleNeutron(eventid);
  //--- Gamma
  //SetParticleGamma1(eventid);
  //--- Electron
  //SetParticleElectron(eventid);
  //--- Cosmic
  //SetParticleCosmic();
  //    if(i%10000==0)    G4cout<<"ParticleGun:"<<i<<G4endl;
  particleGun->GeneratePrimaryVertex( anEvent );
  }

  // G4float rvp = G4UniformRand(); 
  // if(rvp<=0.5){
  //   SetParticleAlpha(eventid);
  //   particleGun->GeneratePrimaryVertex( anEvent );
  // }
  // else{
  //   SetParticleGamma1(eventid);
  //   particleGun->GeneratePrimaryVertex( anEvent );
  // }

  // SetParticleGamma1(eventid);
  // particleGun->GeneratePrimaryVertex( anEvent );

  //  SetParticleGamma2(eventid);
  //  particleGun->GeneratePrimaryVertex( anEvent );

  // double rn = G4UniformRand();
  // if(rn<0.89){
  //   SetParticleGamma1(eventid);
  //   particleGun->GeneratePrimaryVertex( anEvent );
  // }
  // else{
  //   SetParticleGamma2(eventid);
  //   particleGun->GeneratePrimaryVertex( anEvent );
  // }

  // energy_pri = 0.;
  // energy_pri = particleGun->GetParticleEnergy();

}

void DMXPrimaryGeneratorAction::SetParticleAlpha(G4int ievt){

  //---Ion
  G4String particleName;
  G4ParticleDefinition* particle;
  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();

  //alpha---
  particle = particleTable->FindParticle(particleName="alpha");
  particleGun->SetParticleDefinition(particle);

  //----- Set Energy
  // G4float Rvp = G4UniformRand(); 
  // if(0<=Rvp && Rvp<12.7){
  //   particleGun->SetParticleEnergy(5.443*MeV);
  // }
  // else if(12.7<=Rvp && Rvp<1){
  //   particleGun->SetParticleEnergy(5.486*MeV);
  // }
  particleGun->SetParticleEnergy(5.486*MeV);
  //  particleGun->SetParticleEnergy(1.0*MeV);

  G4double cost=0.;

  bool IsSelect=false;

  while(!(IsSelect)){
    //    cost = 1.-2.*(G4double)G4UniformRand();
    cost = 1.-2.*G4UniformRand();
    if(cost>0){
      G4double tant = tan(acos(cost));
      if(tant<1./4) IsSelect=true;
    }
  }

  //  G4double test = G4UniformRand();
  //  G4float test2 = G4UniformRand();

  //  G4cout << "phi1 = " << std::setprecision(15) << 2*3.141592*(G4double)test << G4endl;
  //  G4cout << "phi2 = " << std::setprecision(15) << 2.*3.141592*(G4double)test << G4endl;

  //  const G4double PI = 3.14159265359;

  //----- Set Momentum
  //  G4double cost = 1.-2.*(double)G4UniformRand();
  G4double sint = sqrt((1.-cost)*(1.+cost));
  //  G4double phi  = 2*3.141592*(G4double)G4UniformRand();
  //  G4double phi  = 2*3.141592*G4UniformRand();
  //  G4double phi  = 2*3.1*(G4double)G4UniformRand();
  G4double phi  = 2*PI*G4UniformRand();
  G4double sinp = sin(phi);
  G4double cosp = cos(phi);

  G4float xx=sint*cosp;
  G4float yy=sint*sinp;
  G4float zz=cost;

  // G4double xx=sint*cosp;
  // G4double yy=sint*sinp;
  // G4double zz=cost;

  //  zz = sqrt(pow(zz,2));
  G4ThreeVector vecgun(xx,yy,zz);
  //  G4ThreeVector vecgun(0,0,1);
  particleGun->SetParticleMomentumDirection( vecgun );
  G4float xxx = 1.-2.*G4UniformRand();
  G4float yyy = 1.-2.*G4UniformRand();
  G4float zzz = 1.-2.*G4UniformRand();
  //  G4ThreeVector vecpos(2.*xxx, 2.*yyy, -25.3); // GasSetup Alpha
  G4ThreeVector vecpos(0., 0., -25.3); // GasSetup Alpha
  //  G4ThreeVector vecpos(-5.0, 0., -25.3); // GasSetup Alpha
  //  G4ThreeVector vecpos(0., 0., -27.7); // Alpha Run
  //  G4ThreeVector vecpos(0., 0., 0.); // Alpha Run
  particleGun->SetParticlePosition( vecpos );

  // if(0<=Rvp && Rvp<12.9/(12.9+49.5)){
  //   zz = -sqrt(pow(zz,2));
  //   G4ThreeVector vecgun(xx,yy,zz);
  //   particleGun->SetParticleMomentumDirection( vecgun );
  //   G4ThreeVector vecpos(0., 0., 27.7); // Alpha Run
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // else if(12.9/(12.9+49.5)<=Rvp && Rvp<1){
  //   zz = sqrt(pow(zz,2));
  //   G4ThreeVector vecgun(xx,yy,zz);
  //   particleGun->SetParticleMomentumDirection( vecgun );
  //   G4ThreeVector vecpos(0., 0., -27.7); // Alpha Run
  //   particleGun->SetParticlePosition( vecpos );
  // }
  
  //  G4ThreeVector vecgun(-1.,0.,0.);
  //  G4ThreeVector vecgun(0.,0.,1.);

  // G4float radius = 27.8*mm; 
  // G4float xxx = 1.-2.*G4UniformRand();
  // G4float yyy = 1.-2.*G4UniformRand();
  // G4float zzz = 1.-2.*G4UniformRand();
  //  G4ThreeVector vecpos(0., 0., -27.7); // RunV Alpha
  //  G4ThreeVector vecpos(2.*xxx, 2.*yyy, -25.3); // GasSetup Alpha

}

void DMXPrimaryGeneratorAction::SetParticleIon(G4int ievt){

  //---Ion
  //  G4String particleName;
  //  G4ParticleDefinition* particle;
  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();

  //alpha---
  //particle = particleTable->FindParticle(particleName="alpha");
  //particleGun->SetParticleDefinition(particle);

  //Ar ion---
  G4int iZ = 18;
  G4int iA = 40;
  G4float iC = 1.*eplus;
  G4float excitE = 0.*keV;
  particle = particleTable->GetIon(iZ,iA,excitE);
  particleGun->SetParticleDefinition(particle);
  particleGun->SetParticleCharge(iC);
  
  //----- Set Energy
  particleGun->SetParticleEnergy(300.*keV);

  //----- Set Momentum
  G4double cost = 1.-2.*G4UniformRand();
  G4double sint = sqrt((1.-cost)*(1.+cost));
  //  G4double phi  = 2*3.141592*G4UniformRand();
  G4double phi  = 2*PI*G4UniformRand();
  G4double sinp = sin(phi);
  G4double cosp = cos(phi);

  G4float xx=sint*cosp;
  G4float yy=sint*sinp;
  G4float zz=cost;
  G4ThreeVector vecgun(xx,yy,zz);
  //G4ThreeVector vecgun(0.,0.,-1.);
  particleGun->SetParticleMomentumDirection( vecgun );

  ievt = ievt / 100;
  G4int iz = ievt / 21 % 18;
  G4int ir = ievt % 21;

  G4ThreeVector vecpos(-120.+ir*10.*mm, 0., -545.+iz*10.*mm);
  //  G4ThreeVector vecpos(0., 0., -505.*mm);
  particleGun->SetParticlePosition( vecpos );

  //G4cout<<"Particleposition("<<-120.+ir*10.*mm<<","<<-545.+iz*10.*mm<<")"<<G4endl;

}

void DMXPrimaryGeneratorAction::SetParticlePhoton(G4int ievt,G4int iphoton){

  //---photon
  G4String particleName;
  G4ParticleDefinition* particle;
  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
  
  //---Set Photon Energy
  //  G4float PhotonEnergy = 9.69*eV;
  G4float PhotonEnergy = 9.686*eV;
  //  G4float PhotonEnergy = 1.5*eV;
  //  G4float PhotonEnergy = 2.952*eV;
  //  G4float PhotonEnergy = 3.0996*eV;
  //  PhotonEnergy = GetPhotonEnergy();

  particle = particleTable->FindParticle(particleName="opticalphoton");
  particleGun->SetParticleDefinition(particle);
  particleGun->SetParticleEnergy(PhotonEnergy);

  G4double cost = 1.-2.*G4UniformRand();
  G4double sint = sqrt((1.-cost)*(1.+cost));

  //  G4double phi  = 2*3.141592*G4UniformRand();
  G4double phi  = 2*PI*G4UniformRand();
  G4double sinp = sin(phi);
  G4double cosp = cos(phi);

  G4ThreeVector vecgun(sint*cosp,sint*sinp,cost);
  //  G4ThreeVector vecgun(1,0,1);
  particleGun->SetParticleMomentumDirection( vecgun );

  //**********TPB Kiz***************//

  G4float xxx = 0.01*(1.-2.*G4UniformRand());
  G4float yyy = 0.01*(1.-2.*G4UniformRand());

  //  G4ThreeVector vecpos( rr*cos(pphi), rr*sin(pphi), zz);
  //  G4ThreeVector vecpos( xx, yy, zz);
  //  G4ThreeVector vecpos( 0, -84.12, -64.2); // PreRunIV
  //  G4ThreeVector vecpos(0., 0., -27.68); // PreRunIV
  G4ThreeVector vecpos(0., 0., 0.); // PreRunIV
  //  G4ThreeVector vecpos(0., 0., -25.3); // PreRunIV
  //  G4ThreeVector vecpos(0., 0., -27.7); // PreRunIV
  //  G4ThreeVector vecpos(0., 0., 0.7); // PreRunIV
  //  G4ThreeVector vecpos( xxx, yyy, 0);
  //  G4ThreeVector vecpos( 130., 0, 0);

  //*************Position Scan******************//

  //   G4cout << "iphoton = " << iphoton << G4endl;

  // if(iphoton==0){
  //   G4ThreeVector vecpos(2.5, 0., 27.8/12);
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==1){
  //   G4ThreeVector vecpos(7.5, 0., 27.8/12);
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==2){
  //   G4ThreeVector vecpos(12.5, 0., 27.8/12);
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==3){
  //   G4ThreeVector vecpos(17.5, 0., 27.8/12);
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==4){
  //   G4ThreeVector vecpos(22.5, 0., 27.8/12);
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==5){
  //   G4ThreeVector vecpos(27.5, 0., 27.8/12);
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==6){
  //   G4ThreeVector vecpos(2.5, 0., (27.8/12)+((27.8/6)*1.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==7){
  //   G4ThreeVector vecpos(7.5, 0., (27.8/12)+((27.8/6)*1.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==8){
  //   G4ThreeVector vecpos(12.5, 0., (27.8/12)+((27.8/6)*1.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==9){
  //   G4ThreeVector vecpos(17.5, 0., (27.8/12)+((27.8/6)*1.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==10){
  //   G4ThreeVector vecpos(22.5, 0., (27.8/12)+((27.8/6)*1.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==11){
  //   G4ThreeVector vecpos(27.5, 0., (27.8/12)+((27.8/6)*1.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==12){
  //   G4ThreeVector vecpos(2.5, 0., (27.8/12)+((27.8/6)*2.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==13){
  //   G4ThreeVector vecpos(7.5, 0., (27.8/12)+((27.8/6)*2.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==14){
  //   G4ThreeVector vecpos(12.5, 0., (27.8/12)+((27.8/6)*2.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==15){
  //   G4ThreeVector vecpos(17.5, 0., (27.8/12)+((27.8/6)*2.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==16){
  //   G4ThreeVector vecpos(22.5, 0., (27.8/12)+((27.8/6)*2.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==17){
  //   G4ThreeVector vecpos(27.5, 0., (27.8/12)+((27.8/6)*2.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==18){
  //   G4ThreeVector vecpos(2.5, 0., (27.8/12)+((27.8/6)*3.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==19){
  //   G4ThreeVector vecpos(7.5, 0., (27.8/12)+((27.8/6)*3.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==20){
  //   G4ThreeVector vecpos(12.5, 0., (27.8/12)+((27.8/6)*3.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==21){
  //   G4ThreeVector vecpos(17.5, 0., (27.8/12)+((27.8/6)*3.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==22){
  //   G4ThreeVector vecpos(22.5, 0., (27.8/12)+((27.8/6)*3.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==23){
  //   G4ThreeVector vecpos(27.5, 0., (27.8/12)+((27.8/6)*3.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==24){
  //   G4ThreeVector vecpos(2.5, 0., (27.8/12)+((27.8/6)*4.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==25){
  //   G4ThreeVector vecpos(7.5, 0., (27.8/12)+((27.8/6)*4.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==26){
  //   G4ThreeVector vecpos(12.5, 0., (27.8/12)+((27.8/6)*4.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==27){
  //   G4ThreeVector vecpos(17.5, 0., (27.8/12)+((27.8/6)*4.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==28){
  //   G4ThreeVector vecpos(22.5, 0., (27.8/12)+((27.8/6)*4.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==29){
  //   G4ThreeVector vecpos(27.5, 0., (27.8/12)+((27.8/6)*4.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==30){
  //   G4ThreeVector vecpos(2.5, 0., (27.8/12)+((27.8/6)*5.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==31){
  //   G4ThreeVector vecpos(7.5, 0., (27.8/12)+((27.8/6)*5.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==32){
  //   G4ThreeVector vecpos(12.5, 0., (27.8/12)+((27.8/6)*5.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==33){
  //   G4ThreeVector vecpos(17.5, 0., (27.8/12)+((27.8/6)*5.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==34){
  //   G4ThreeVector vecpos(22.5, 0., (27.8/12)+((27.8/6)*5.));
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // if(iphoton==35){
  //   G4ThreeVector vecpos(27.5, 0., (27.8/12)+((27.8/6)*5.));
  //   particleGun->SetParticlePosition( vecpos );
  // }

  //********************************************//

  particleGun->SetParticlePosition( vecpos );

}

void DMXPrimaryGeneratorAction::SetParticleNeutron(G4int ievt){

  //---Neutron
  G4String particleName;
  G4ParticleDefinition* particle;
  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();

  particle = particleTable->FindParticle(particleName="neutron");
  particleGun->SetParticleDefinition(particle);

  //---Set Energy
  G4float NeutronEnergy;
  //NeutronEnergy = Get241AmBeEnergy()*MeV; // 241Am-Be(Bare)
  //NeutronEnergy = Get252CfEnergy()*MeV;   // 252Cf

  NeutronEnergy = 1.*MeV;
  particleGun->SetParticleEnergy(NeutronEnergy);

  //---Set Momentum
  G4double cost = 1.-2.*G4UniformRand();
  G4double sint = sqrt((1.-cost)*(1.+cost));
  G4double phi  = 2*PI*G4UniformRand();
  //  G4double phi  = 2*3.141592*G4UniformRand();
  G4double sinp = sin(phi);
  G4double cosp = cos(phi);

  //  G4ThreeVector vecgun(sint*cosp,sint*sinp,cost);
  G4ThreeVector vecgun(1.,0.,0.);
  particleGun->SetParticleMomentumDirection( vecgun );

  //---Set Position
  G4int iz = ievt / 1000;
  //  G4float xx = -172.*mm;
  G4float xx = -180.*mm;
  G4float yy = 0.*mm;
  G4float zz = (-545.+iz*20.)*mm;

  //G4ThreeVector vecpos(xx, yy, zz);
  //  G4cout<<"Particleposition("<<xx<<","<<yy<<","<<zz<<")"<<G4endl;
  G4cout<<"NeutronEnergy = "<<G4BestUnit(NeutronEnergy,"Energy")<<G4endl;
  //  G4ThreeVector vecpos(xx, yy, -455.*mm);
  G4ThreeVector vecpos(0., 0., 0.);

  particleGun->SetParticlePosition( vecpos );
}

void DMXPrimaryGeneratorAction::SetParticleGamma1(G4int ievt){
  //---Gamma
  G4String particleName;
  G4ParticleDefinition* particle;
  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();

  particle = particleTable->FindParticle(particleName="gamma");
  particleGun->SetParticleDefinition(particle);

  //---Set Energy
  G4float GammaEnergy;
  //  GammaEnergy = 1.33*MeV;
  //  GammaEnergy = 0.622*MeV;
  //  GammaEnergy = 0.511*MeV;
  GammaEnergy = 0.122*MeV;
  //  GammaEnergy = 0.059*MeV;
  particleGun->SetParticleEnergy(GammaEnergy);

  //---Set Momentum
  G4double cost = 1.-2.*G4UniformRand();
  G4double sint = sqrt((1.-cost)*(1.+cost));
  //  G4double phi  = 2*3.14159265358979*G4UniformRand();
  G4double phi  = 2*PI*G4UniformRand();
  //  G4double phi  = 0.5*3.141592*(1.-2.*G4UniformRand());
  G4double sinp = sin(phi);
  G4double cosp = cos(phi);

  G4float xx=sint*cosp;
  G4float yy=sint*sinp;
  G4float zz=cost;

  // G4float RVP = G4UniformRand(); 

  // if(RVP<0.5){
  //   zz = -sqrt(pow(zz,2));
  //   G4ThreeVector vecgun(xx,yy,zz);
  //   particleGun->SetParticleMomentumDirection( vecgun );
  //   G4ThreeVector vecpos(0., 0., 27.7); // Alpha Run
  //   particleGun->SetParticlePosition( vecpos );
  // }
  // else{
  //   zz = sqrt(pow(zz,2));
  //   G4ThreeVector vecgun(xx,yy,zz);
  //   particleGun->SetParticleMomentumDirection( vecgun );
  //   G4ThreeVector vecpos(0., 0., -27.7); // Alpha Run
  //   particleGun->SetParticlePosition( vecpos );
  // }

  yy = sqrt(pow(yy,2));
  //  G4ThreeVector vecgun(sint*cosp,sint*sinp,cost);
  G4ThreeVector vecgun(xx,yy,zz);
  //  G4ThreeVector vecgun(0.,1.,0.);
  particleGun->SetParticleMomentumDirection( vecgun );

  // //---Set Position
  // G4int iz = ievt / 10000;
  // G4float PHI = -67.5*deg;
  // G4float xx = 180.*cos(PHI)*mm;
  // G4float yy = 180.*sin(PHI)*mm;
  // G4float zz = -465.*mm;
  // //G4float zz = (-545.+iz*20.)*mm;

  // //G4ThreeVector vecpos(xx, yy, zz);
  // //  G4cout<<"Particle Position("<<xx<<","<<yy<<","<<zz<<")"<<G4endl;
  // //  G4cout<<"      Beam Vector("<<sint*cosp<<","<<sint*sinp<<","<<cost<<")"<<G4endl;
  // //  G4cout<<"Gamma Energy = "<<G4BestUnit(GammaEnergy,"Energy")<<G4endl;
  // //  G4ThreeVector vecpos(xx, yy, zz);
  //  G4ThreeVector vecpos(0.,-165., -64.2);
  G4ThreeVector vecpos(0.,-165., 0.);
  // //  G4ThreeVector vecpos(165.+100., 0., 0.);
  // //  G4ThreeVector vecpos(0., 165., 0.);
  // //  G4ThreeVector vecpos(0., 0., 0.);

  particleGun->SetParticlePosition( vecpos );

}

void DMXPrimaryGeneratorAction::SetParticleGamma2(G4int ievt){
  //---Gamma
  G4String particleName;
  G4ParticleDefinition* particle;
  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();

  particle = particleTable->FindParticle(particleName="gamma");
  particleGun->SetParticleDefinition(particle);

  //---Set Energy
  G4float GammaEnergy;
  //  GammaEnergy = 1.17*MeV;
  GammaEnergy = 0.136*MeV;
  particleGun->SetParticleEnergy(GammaEnergy);

  //---Set Momentum
  G4double cost = 1.-2.*G4UniformRand();
  G4double sint = sqrt((1.-cost)*(1.+cost));
  //  G4double phi  = 2*3.14159265358979*G4UniformRand();
  G4double phi  = 2*PI*G4UniformRand();
  //  G4double phi  = 0.5*3.141592*(1.-2.*G4UniformRand());
  G4double sinp = sin(phi);
  G4double cosp = cos(phi);

  G4float xx=sint*cosp;
  G4float yy=sint*sinp;
  G4float zz=cost;

  yy = sqrt(pow(yy,2));
  G4ThreeVector vecgun(xx,yy,zz);
  //  G4ThreeVector vecgun(sint*cosp,sint*sinp,cost);
  //  G4ThreeVector vecgun(-1.,0.,0.);
  particleGun->SetParticleMomentumDirection( vecgun );

  // //---Set Position
  // G4int iz = ievt / 10000;
  // G4float PHI = -67.5*deg;
  // G4float xx = 180.*cos(PHI)*mm;
  // G4float yy = 180.*sin(PHI)*mm;
  // G4float zz = -465.*mm;
  // //G4float zz = (-545.+iz*20.)*mm;

  //G4ThreeVector vecpos(xx, yy, zz);
  //  G4cout<<"Particle Position("<<xx<<","<<yy<<","<<zz<<")"<<G4endl;
  //  G4cout<<"      Beam Vector("<<sint*cosp<<","<<sint*sinp<<","<<cost<<")"<<G4endl;
  //  G4cout<<"Gamma Energy = "<<G4BestUnit(GammaEnergy,"Energy")<<G4endl;
  //  G4ThreeVector vecpos(xx, yy, zz);
  //  G4ThreeVector vecpos(165., 0., 0.);
  //  G4ThreeVector vecpos(0.,-165., -64.2);
  G4ThreeVector vecpos(0.,-165.,0.);
  //  G4ThreeVector vecpos(0., 165., 0.);
  //  G4ThreeVector vecpos(0., 0., 0.);

  particleGun->SetParticlePosition( vecpos );

}

void DMXPrimaryGeneratorAction::SetParticleElectron(G4int ievt){
  //---Electron
  G4String particleName;
  G4ParticleDefinition* particle;
  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();

  particle = particleTable->FindParticle(particleName="e-");
  particleGun->SetParticleDefinition(particle);

  //---Set Energy
  G4float ElectronEnergy;
  ElectronEnergy = 5.*MeV;
  particleGun->SetParticleEnergy(ElectronEnergy);

  //---Set Momentum
  G4double cost = 1.-2.*G4UniformRand();
  G4double sint = sqrt((1.-cost)*(1.+cost));
  //  G4double phi  = 2*3.14159265358979*G4UniformRand();
  G4double phi  = 2*PI*G4UniformRand();
  //  G4double phi  = 0.5*3.141592*(1.-2.*G4UniformRand());
  G4double sinp = sin(phi);
  G4double cosp = cos(phi);

  G4ThreeVector vecgun(sint*cosp,sint*sinp,cost);
  //  G4ThreeVector vecgun(0.,1.,0.);
  particleGun->SetParticleMomentumDirection( vecgun );

  //---Set Position
  G4int iz = ievt / 10000;
  G4float PHI = -67.5*deg;
  G4float xx = 180.*cos(PHI)*mm;
  G4float yy = 180.*sin(PHI)*mm;
  G4float zz = -465.*mm;
  //G4float zz = (-545.+iz*20.)*mm;

  //G4ThreeVector vecpos(xx, yy, zz);
  //  G4cout<<"Particle Position("<<xx<<","<<yy<<","<<zz<<")"<<G4endl;
  //  G4cout<<"      Beam Vector("<<sint*cosp<<","<<sint*sinp<<","<<cost<<")"<<G4endl;
  //  G4cout<<"Electron Energy = "<<G4BestUnit(ElectronEnergy,"Energy")<<G4endl;
  //  G4ThreeVector vecpos(xx, yy, zz);
  //  G4ThreeVector vecpos(0.,-165., -64.2);
  //  G4ThreeVector vecpos(165.+100., 0., 0.);
  //  G4ThreeVector vecpos(0., 165., 0.);
  G4ThreeVector vecpos(0., 0., 0.);

  particleGun->SetParticlePosition( vecpos );

}

void DMXPrimaryGeneratorAction::SetParticleCosmic(){

  //  gRandom->SetSeed(time(NULL));
  // TFile *fread=new TFile("./inputFiles/MuonFlux.root");
  // TH1F *hSpectrumMu_p=(TH1F*)fread->Get("hSpectrumMu_p");
  // TH1F *hSpectrumMu_n=(TH1F*)fread->Get("hSpectrumMu_n");
  // TH1F *hSpectrumMu  =(TH1F*)fread->Get("hSpectrumMu");

  //  RatioPtoALL = hSpectrumMu_p->Integral()/hSpectrumMu->Integral();

  //---Cosmic

  G4String particleName;
  G4ParticleDefinition* particle;
  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();

  //  float vRandom = gRandom->Rndm();
  //  float MuonEnergy = 0;

  float MuonEnergy = 800.;

  particle = particleTable->FindParticle(particleName="mu+");

  //  if( RatioPtoALL>vRandom ){
  //    particle = particleTable->FindParticle(particleName="mu+");
  //    MuonEnergy = hSpectrumMu_p->GetRandom();
  //  }
  //  else{
  //    particle = particleTable->FindParticle(particleName="mu-");
  //    MuonEnergy = hSpectrumMu_n->GetRandom();
  //  }

  particleGun->SetParticleDefinition(particle);
  particleGun->SetParticleEnergy(MuonEnergy*MeV);

  //  G4cout<<"Ratio = "<<RatioPtoALL<<G4endl;
  //  G4cout<<"MuonType   = "<<particleName<<G4endl;
  //  G4cout<<"MuonEnergy = "<<G4BestUnit(MuonEnergy,"Energy")<<G4endl;

  //  G4float phi  = 2*3.14159265358979*G4UniformRand();
  G4float phi  = 2*PI*G4UniformRand();
  G4float sinp = sin(phi);
  G4float cosp = cos(phi);

  G4float r = 50.*G4UniformRand();

  G4float x = r*cosp;
  G4float y = r*sinp;
  G4float z = 50.*(1.-2.*G4UniformRand());

  // G4float pphi  = 2*3.14159265358979*G4UniformRand();
  // G4float ssinp = sin(pphi);
  // G4float ccosp = cos(pphi);

  // G4float rr = 50.*G4UniformRand();

  // G4float xx = rr*ccosp;
  // G4float yy = rr*ssinp;
  // G4float zz = 50.*(1.-2.*G4UniformRand());

  // G4float xxx;
  // G4float yyy;
  // G4float zzz;

  // if(zz-z<0){
  //   xxx = xx-x;
  //   yyy = yy-y;
  //   zzz = zz-z;
  // }
  // else if(zz-z>=0){
  //   xxx = x-xx;
  //   yyy = y-yy;
  //   zzz = z-zz;
  // }

  TF1 *f = new TF1("f","pow(cos(x),2)",-3.14159265358979,3.14159265358979);

  G4double value=-999;
  G4double valuel=999;
  G4double theta=0.;
  G4double cost=0.;

  while(valuel>value){
    cost = 1.-2.*G4UniformRand();
    theta = acos(cost);
    value = f->Eval(theta);
    valuel = G4UniformRand();
  }

  G4double ccost = -sqrt(pow(cost,2));

  //---Set Momentum
  //  G4double ccost = 1.-2.*G4UniformRand();
  //  G4double ccost = -1.*G4UniformRand();
  G4double ssint = sqrt((1.-ccost)*(1.+ccost));
  //  G4double pphi  = 2*3.14159265358979*G4UniformRand();
  G4double pphi  = 2*PI*G4UniformRand();
  G4double ssinp = sin(pphi);
  G4double ccosp = cos(pphi);

  //  G4ThreeVector vecgun(sint*cosp,sint*sinp,cost);

  G4float xxx = ssint*ccosp;
  G4float yyy = ssint*ssinp;
  G4float zzz = ccost;

  G4float A = pow(xxx,2)+pow(yyy,2)+pow(zzz,2);
  G4float B = x*xxx+y*yyy+z*zzz;
  G4float C = pow(x,2)+pow(y,2)+pow(z,2)-pow(500.,2);

  G4float t = (-B-sqrt(pow(B,2)-A*C))/A;

  G4float sx = x+(t*xxx);
  G4float sy = y+(t*yyy);
  G4float sz = z+(t*zzz);

  G4ThreeVector vecgun(xxx,yyy,zzz);
  //  G4ThreeVector vecgun(-1,0,-1);
  particleGun->SetParticleMomentumDirection( vecgun );

  G4ThreeVector vecpos(sx, sy, sz);
  //  G4ThreeVector vecpos(500, 0, 500);
  //  G4ThreeVector vecpos(0, 0, 500.);
  particleGun->SetParticlePosition( vecpos );

//   fread->Delete();
//   hSpectrumMu_p->Delete();
//   hSpectrumMu_n->Delete();
//   hSpectrumMu  ->Delete();

}

G4double DMXPrimaryGeneratorAction::GetPhotonEnergy(){
  TF1 *fAr=new TF1("fAr","exp(-0.5*((x-128.0)/2.929)*((x-128.0)/2.929))",0,400.);
  G4double PlankH = 4.135667516e-15*eV*s;
  G4double LightC = 2.99792458e8*m/s; 
  G4double PhotonEnergy = PlankH*LightC/(fAr->GetRandom()*nm);
  fAr->Delete();

  return PhotonEnergy;
}
G4double DMXPrimaryGeneratorAction::Get252CfEnergy(){
  TF1 *fNspectrum=new TF1("fNspectrum","sqrt(x)*exp(-x/1.42)",0.01,100);
  G4double NeutronEnergy = fNspectrum->GetRandom();
  fNspectrum->Delete();

  return NeutronEnergy;
}
G4double DMXPrimaryGeneratorAction::Get241AmBeEnergy(){
  TH1F *h1=new TH1F("h1","h1",10000,0,100);

  ifstream fin("./inputFiles/AmBeTable.dat");
  if(!fin){G4cout<<"No Such File"<<G4endl;return 0;}

  while( !fin.eof() ){
    float Energy,Cf,AmBe,tmp;
    fin >> Energy >> Cf >> tmp >> AmBe >> tmp;
    h1->Fill((Energy/1.0e6),AmBe);
  }
  int nbin = h1->GetNbinsX();
  for(int ibin=0;ibin<nbin;ibin++){
    float vbin = h1->GetBinContent(ibin+1);
    if( ibin<15.8/(100/10000.)&&!vbin&&ibin>0 ){
      h1->SetBinContent(ibin+1,h1->GetBinContent(ibin));
    }
  }
  G4double NeutronEnergy = h1->GetRandom();
  h1->Delete();

  return NeutronEnergy;
}
